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ABSTRACT 

The  Kirchhoff  integral  provides  a  robust  method  for  implementing  seismic 
modeling  and  prestack  depth  migration,  which  can  handle  lateral  velocity  varia¬ 
tion  and  turning  waves.  With  a  little  extra  computation  cost,  the  Kirchhoff-type 
migration  can  obtain  multiple  outputs  that  have  the  same  phase  but  different 
amplitudes,  compared  with  that  of  other  migration  methods.  The  ratio  of  these 
amplitudes  is  helpful  in  computing  some  quantities  such  as  reflection  angle. 

Here,  I  develop  a  seismic  modeling  and  prestack  depth  migration  method 
based  on  the  Kirchhoff  integral,  that  handles  both  laterally  variant  velocity  and 
a  dip  beyond  90  degrees.  The  method  uses  a  finite-difference  algorithm  to  cal¬ 
culate  traveltimes  and  WKBJ  amplitudes  for  the  Kirchhoff  integral.  Compared 
to  ray-tracing  algorithms,  the  finite-difference  algorithm  gives  an  efficient  imple¬ 
mentation  and  single-valued  quantities  (first  arrivals)  on  output.  In  my  finite- 
difference  algorithm,  the  upwind  scheme  is  used  to  calculate  traveltimes,  and  the 
Crank-Nicolson  scheme  is  used  to  calculate  amplitudes.  Moreover,  interpolation 
is  applied  to  save  computation  cost. 

The  modeling  and  migration  algorithms  here  require  a  smooth  velocity  func¬ 
tion.  I  develop  a  velocity-smoothing  technique  based  on  damped  least-squares  to 
aid  in  obtaining  a  successful  migration.  This  velocity-smoothing  technique  also 
can  be  used  to  improve  results  of  other  migration  algorithms,  such  as  Gaussian 
beam  migration. 


INTRODUCTION 

Seismic  modeling  and  migration  play  an  important  role  in  seismic  data  processing. 
To  treat  complex  media,  one  needs  prestack  depth  migration  that  can  handle  lateral 
velocity  variation  and  reflector  dips.  Conventional  techniques,  such  as  the  down¬ 
ward  continuation  of  sources  and  geophones  by  finite  difference  (S-G  finite-difference 
migration),  are  relatively  slow  and  dip-limited.  Compared  to  S-G  finite-difference 
migration,  the  Kirchhoff  integral  implements  prestack  migration  relatively  efficiently 
and  has  no  dip  limitation. 
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Bleistein  et  al.  (1987)  derived  the  Kirchhoff  integral  method  by  using  the  WKBJ 
approximation.  This  method  treats  amplitude  in  migration  in  a  WKBJ-consistent 
manner  so  that  the  output  is  the  reflectivity  function.  Furthermore.  Bleistein  showed 
that  one  can  design  weights  in  the  Kirchhoff  integral  to  determine  quantities  such 
as  the  reflection  angle,  based  on  the  stationary-phase  principle.  In  Bleistein’s  inte¬ 
gral  formula,  the  integrand  consists  of  traveltime,  WKBJ  amplitude,  ray  parameters, 
and  so  on.  Two  commonly  used  approaches  to  calculate  traveltimes  and  the  other 
quantities  in  the  Kirchhoff  integral  are  rav-tracing  and  finite-difference  applied  to 
the  eikonal  equation.  The  ray-tracing  algorithm  calculates  traveltimes  and  ampli¬ 
tudes  on  each  ray.  Usually,  the  rays  do  not  pass  through  output  grid  points,  so  we 
need  to  interpolate  to  obtain  traveltimes  and  amplitudes  at  the  grid  locations.  This 
interpolation  is  complicated  when  the  raypaths  have  a  caustic,  and  thus  limits  the 
efficiency  of  the  ray-tracing  algorithm.  In  contrast,  the  finite-difference  algorithm  cal¬ 
culates  traveltimes  and  amplitudes  directly  at  output  grid  locations.  For  caustics,  the 
finite-difference  algorithm  calculates  the  first-arrival  traveltime  and  the  corresponding 
amplitudes. 

Traveltimes  satisfy  the  eikonal  equation,  and  amplitude  terms  satisfy  linear  par¬ 
tial  differential  equations  that  depend  on  traveltime  derivatives.  Van  Trier  and  Svmes 
(1990)  introduced  a  finite-difference  scheme  for  solving  this  equation  to  obtain  the 
traveltimes  and  traveltime  derivatives.  Pusey  and  Vidale  (1991)  used  a  similar,  ex¬ 
plicit  scheme  to  solve  for  the  WKBJ  amplitudes.  For  the  explicit  scheme,  one  has  to 
choose  a  small  enough  step  size  for  a  stability  when  the  velocity  function  is  varied.  In 
this  paper,  I  use  the  Crank-Nicolson  scheme  to  solve  for  the  amplitudes.  Compared 
to  the  explicit  scheme,  the  Crank-Nicolson  scheme  is  second-order  accurate  and  ab¬ 
solutely  stable  so  that  computation  cost  is  made  relatively  small  for  variable  velocity 
by  choosing  large  step  sizes. 

Calculation  cost  of  pointwise  traveltimes  and  amplitudes  by  the  above  method 
would  dominant  the  total  calculation  cost  of  the  Kirchhoff  integral  method.  For  each 
source  or  receiver,  the  calculation  cost  depends  on  the  number  of  grid  points.  To 
speed  up  the  Kirchhoff  integral  method,  I  apply  interpolation  to  avoid  computation 
on  the  entire  grid. 

A  smooth  velocity  function  is  required  by  most  migration  algorithms.  A  few 
methods  for  smoothing  the  model  velocity  have  been  developed  in  the  past.  The 
windowed-averaging  method,  for  example,  calculates  a  smooth  velocity  value  at  one 
point  by  averaging  velocity  values  over  a  depth-midpoint  window  centered  at  this 
point.  If  the  original  velocity  is  discontinuous,  it  is  easy  to  show  that  the  smoothed 
velocity  from  the  windowed-averaging  method  does  not  have  the  continuity  of  the 
first  derivative.  Consequently,  this  method  is  not  effective  when  the  velocity  has  a 
strong  discontinuity.  Here,  I  present  a  velocity-smoothing  technique  based  on  damped 
least-squares.  The  curvature  of  the  raypath  is  a  key  factor  in  migration  imaging.  The 
smaller  the  curvature,  the  more  stable  the  migration  result.  I  derive  a  representation 
for  the  curvature  that  depends  on  first  derivatives  of  velocity.  A  smooth  velocity 
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function  is  sought  that  minimizes  the  weighted  sum  of  (1)  the  deviation  between  the 
smooth  velocity  and  the  original  one.  and  (2)  the  first  derivatives  of  velocity  with 
respect  to  spatial  variables.  I  use  the  result  that  the  reciprocal  of  the  curvature  of 
the  raypath  is  proportional  to  the  velocity  gradient.  Therefore,  my  method  may 
remove  migration  artifacts  by  suppressing  the  curvature  of  the  raypath.  Moreover, 
this  method  allows  local  variation  in  the  degree  of  smoothing  allowed  by  using  a 
window  function,  so  that  the  velocity  is  smoothed  more  in  rougher  areas. 


FINITE-DIFFERENCE  ALGORITHM 

The  Kirchhoff  integral  method  can  be  represented  by 


output  =  integral{weight  ■  input}. 


For  modeling  problems,  the  input  consists  of  velocity  layers,  and  the  output  is  seismic 
traces;  for  inversion  problems,  the  input  is  seismic  traces,  the  output  is  a  structural 
image.  Both  require  common  quantities  to  be  calculated.  These  quantities  include 
(Bleistein,  1986): 


r 

© 

o 

P 

dd/dx 


traveltime, 

propagation  angle, 

running  ray  parameter, 

incident  angle  from  source  or  reciever, 

geometrical  spreading  parameter. 


For  each  source  or  receiver,  the  above  quantities  satisfy  the  follow  eqations; 
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Equation  (1)  is  the  eikonal  equation.  Equations  (3)  and  (4)  are  derived  by  Pusey 
and  Vidale  (1991).  Equation  (5)  follows  from  equation  (3).  Van  Trier  et  al.  (1990) 
used  the  efficient  upwind  scheme  to  solve  this  equation  for  traveltime  and  traveltime 
derivatives.  However,  computation  in  Cartesian  coordinates  requires  that  dr/dz  be 
positive,  so  turning  waves  cannot  be  handled.  Computation  in  polar  coordinates  does 
not  suffer  this  limitation.  Let  us  introduce  the  coordinate  transform 


r  =  x,  +  r  sin  0,  z  =  rcos9. 


In  polar-coordinates,  equations  (1)  through  (5)  become 


/ dr\2  1  / dr\ 2 _  1 

\dry  v2(r,9)' 


(6) 


where 
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Although  more  computation  is  required  in  polar-coordinates,  the  boundary  conditions 
are  more  easily  treated  and,  more  importantly,  we  can  handle  turning  waves. 


The  coefficient  p  of  equation  (10)  is  inside  the  differential  operator.  If  an  error 
exists  in  numerical  computation  of  p,  the  solution  of  equation  (10)  may  suffer.  Instead 
of  solving  equation  (10),  I  use  an  approximate  formula 


dd  v0r 

—  =  — COS  (0-0),  (11) 

where  vq  is  the  velocity  value  at  the  source  position.  Formula  (11),  derived  in 
Appendix  A,  is  exact  for  linear  velocity.  When  velocity  is  strongly  nonlinear,  the 
geometry-spreading  factor  based  on  formula  (11)  is  not  accurate.  However,  not  only 
the  geometry-spreading  factor  but  also  the  finite  difference  algorithm,  which  requires 
that  dr /dr  be  positive,  suffers  from  strong-nonlinear  velocity  .  The  smoothing  tech¬ 
nique  in  this  paper  will  produce  a  sufficiently  smooth  velocity  to  solve  these  two 
problems. 
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Crank-Nicolson  scheme 

Van  Trier  et  al.  (1990)  gives  an  upwind  finite-difference  scheme  for  solving  equa 
tion  (6).  Also  that  scheme  yields  derivative  values  for  use  in  equation  (7).  Pusey  and 
Vidale  (1991)  use  a  similar  scheme  to  solve  equation  (8).  Those  upwind  schemes  are 
explicit.  Therefore,  one  has  to  choose  a  small  enough  step  size  for  stability  when  the 
velocity  function  is  varied.  In  this  paper,  the  Crank-Nicolson  scheme  is  applied  to 
equations  (8)  and  (9).  Compared  to  the  explicit  scheme,  the  Crank  Nicolson  scheme 
is  second-order  accurate  and  absolutely  stable  so  that  computation  cost  is  made  rel¬ 
atively  small  for  variable  velocity  by  choosing  large  steps. 


Let 


dr 


_  1  dr 

q~7*d0' 


With  the  Crank-Nicolson  approximation,  equation  (8)  becomes 


^+l - W / j+i ,  t  1 _ «j+i 
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and  equation  (9)  becomes 
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(13) 


where  j  is  the  index  of  r,  i  is  the  index  of  0,  and  Ar  and  Ad  define  the  grid  cell  sizes. 

For  each  j,  equations  (12)  and  (13)  are  tridiagonal  systems  that  can  be  solved 
relatively  efficiently. 


INTERPOLATION  TECHNIQUE 

Calculation  cost  of  pointwise  traveltimes  and  amplitudes  by  the  above  method 
would  dominant  the  total  calculation  cost  of  the  KirchhofF  integral  method.  For  each 
source  or  receiver,  the  calculation  cost  depends  on  the  number  of  grid  points.  To 
speed  up  the  KirchhofF  integral  method,  I  apply  interpolation  at  two  stages  so  as  to 
avoid  computation  on  the  entire  grid. 

One  stage  is  the  interpolation  between  the  sources  or  receivers:  I  calculate  trav¬ 
eltimes  and  amplitudes  only  at  selected  sources  or  receivers,  and  then  interpolate 
traveltimes  and  amplitudes  at  the  other  sources  or  receivers.  Suppose  that  the  trav¬ 
eltime  functions  r(x,  z\ x3l )  and  t(x,z;x3j)  from  two  sources  x3,  and  x3j  have  been 
calculated.  I  interpolate  a  traveltime  function  from  source  x,  (x3l  <  x3  <  x32)  by 

t{x,z,x3)  =  Xt{x+  x32-x3,z-,x3j)  +  (1-  A)t(x+£Si  —xa,  2;xSl).,  (14) 

where  A  =  (x3  —  xai)/(x32  —  x4,).  When  the  velocity  is  only  depth-dependent,  this 
linear  interpolation  is  exact  because  the  traveltime  is  invariant  for  parallel  movement; 
i.e.,  for  any  h, 

r(x  +  h ,  2;  x,  +  h)  =  r(x,  2;  x3). 
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For  a  general  velocity  function  v(x,z),  the  interpolation  error  depends  on  dv/dx. 
Reducing  this  error  is  one  motivation  for  smoothing  the  velocity. 

The  other  interpolation  is  between  grid  points.  To  save  computation,  I  use  a 
relatively  coarse  grid  in  the  traveltime  and  amplitude  calculation.  Then,  to  preserve 
resolution,  I  use  a  finer  grid  for  migration  output  with  traveltimes  and  amplitudes 
calculated  by  bilinear  interpolation. 


VELOCITY  SMOOTHING 


A  smooth  velocity  function  is  required  by  most  migration  algorithms.  Before  a 
smoothing  technique  is  selected,  it  is  essential  to  know  how  smoothness  of  velocity 
affects  migration  results. 

Numerous  migration  algorithms  are  based  on  ray  theory  (high-frequency  assump¬ 
tion).  The  smaller  the  curvature  or  the  larger  the  radius  of  curvature,  the  better  is 
the  quality  of  the  raypath.  Generally,  a  large  curvature  increases  the  sensitivity  of 
a  raypath  to  the  incident  angle  from  the  source,  resulting  in  sparse  raypath  cover¬ 
age.  In  addition,  a  large  curvature  may  cause  troubles  in  the  imaging  process.  In 
Gaussian  beam  migration,  for  example,  traveltime  at  one  point  in  the  vicinity  of  a 
central  ray  are  calculated  by  projecting  this  point  onto  the  ray  and  then  by  using 
the  Taylor  series  expansion.  If  the  curvature  of  the  ray  is  too  large,  the  projecting 
point  is  not  unique  and  also  the  Taylor  series  expansion  is  inaccurate,  so  that  a  poor 
migration  result  may  suffer  from  the  inaccurate  calculation  of  the  traveltimes.  For 
another  example,  the  finite-difference  method  in  this  paper  will  fail  if,  at  any  point 
of  the  raypath,  there  is  more  than  a  90-degree  difference  between  the  ray  direction 
and  the  direction  from  the  source  position  to  this  point.  Therefore,  control  of  the 
curvature  of  the  raypath  is  necessary.  From  Dohr  (1985,  p.  23),  the  curvature  k  can 
be  represented  by 

_  1  _  dv  cos  9  dv  sin  0 

R  dx  v  dz  v 

where  R  is  the  radius  of  the  curvature,  v  is  the  velocity,  and  9  is  the  propagation 
angle.  Equation  (15)  shows  that  the  first  derivatives  dominate  the  curvature  of  the 
ray.  Based  on  this  characteristic,  I  used  the  damped  least-squares  method  to  calculate 
a  smooth  velocity  that  has  suppressed  the  first  derivatives  and,  therefore,  reduced  the 
raypath  curvature.  Use  of  a  smooth  velocity  function  is  also  required  in  the  WKBJ 
approximation  and  in  the  interpolation  for  traveltimes. 


by 


Let  v(x,z)  be  the  original  velocity;  then  a  smooth  velocity  vs{x,  z)  is  determined 


(16) 


(dv  \2 

j  ( v3(x ,  z)  —  v(x,  z ))2  dx  +  o?x  J  w(x,  z)  f  ~  J  dx  =  min, 

J K(x,  z)  -  v3(x,  z))2  dz  +  a2  J  w{x,  z)  dz  =  min’  (17) 
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where  ax,az  are  the  smoothing  parameters  and  w(x,z)  is  a  window  function.  The 
window  function,  ranging  from  0  to  1 ,  allows  local  variation  in  the  degree  of  smoothing 
desired.  Roughness  of  the  original  velocity  function  is  usually  not  uniform,  so  one 
can  design  a  window  function  such  that  velocity  is  smoothed  more  in  rougher  areas. 

Equation  (16)  smooths  the  original  velocity  along  the  x-direction,  from  which 
va(x,z)  is  obtained;  equation  (17)  smooths  va(x,z)  along  the  2-direction  to  get  the 
smoothed  velocity  va(x,z).  In  the  wavenumber  domain  and  for  w(x,z)  =  1,  the 
smooth  velocity  can  be  determined  by 


VM(k„k,) 


V(kx,kz) 

(l+Q2fc2)(l+a2fc2)’ 


(18) 


where  Va  and  V  are  the  Fourier  transforms  of  va  and  v  respectively.  Therefore,  large- 
wavenumber  components  are  suppressed.  In  discretization,  equation  (18)  becomes 


V3(kz,k2) 


_ V(kM _ 

(1  +  4(aI/Ax)2sin2(A:;rAx/2))  (1  +  4(az/ Az)2  sin2  (k.Az/ 2))  ’ 


(19) 


where  Ax  and  A z  define  the  grid-cell  sizes.  A  proof  of  formula  (18)  and  formula  (19) 
is  in  Appendix  B.  I  define  the  unit  of  the  smoothing  parameters  by 


or  =  Ax/2,  az  =  Az/2. 

For  such  smothing  parameters,  the  Nyquist-wavenumber  components  in  x  and  2  di¬ 
rections,  respectively,  reduce  by  1/2,  i.e., 

V3(2/Ax,  2/A2)  =  ~V(2/Ax,2/A2). 

The  larger  the  values  of  ax  and  a.,  the  smoother  will  be  v3(x,  2)  and  the  larger  will 
be  the  difference  between  v(x,  2)  and  v,(x,  2).  The  following  formula  may  be  used  to 
measure  the  difference  between  v(x,z)  and  v,(x,  2): 

2/  x  /  Jo  (».(*. z')  ~  v(x’  z'))2dz'dx 
E  - '  ( 0) 

Based  on  my  experience, if  e(z)  is  greater  than  0.1,  the  smoothing  parameters  are  too 
large — the  velocity  is  oversmoothed. 


COMPUTER  IMPLEMENTATION 

I  applied  the  modeling  and  migration  method  to  synthetic  models.  Also,  I  tested 
the  effect  of  velocity  smoothing  on  the  Kirchhoff  migration  of  this  paper  and  Gaussian 
beam  migration. 
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Synthetic  Data 

The  first  model  shown  in  Figure  la,  consists  of  five  reflectors,  each  with  a  dipping 
and  horizontal  segment.  Dips  for  the  dipping  segments  range  from  30  to  90  degrees  in 
15-degree  increments.  The  velocity  function  consists  of  a  linear  function.  r( z)  =  1.5+2 
km/s,  plus  a  lateral  variation  that  is  cos(^(x  —  1.5))  multiplied  by  a  depth-dependent 
function,  as  shown  in  Figure  lb.  The  maximum  lateral  variation  is  12  percent  at 
r  =  1  km.  With  the  model  and  the  velocity  function,  I  generated  two  data  sets: 
zero-offset  and  an  offset  of  1000  m.  Figures  2  and  3  show  the  modeling  data  and 
inversion  results  for  the  two  offsets.  Because  of  the  lateral  variation  in  velocity,  even 
horizontal  reflectors  give  curve  events  in  the  data.  After  migration,  all  reflectors, 
horizontal  through  vertical,  are  imaged  to  their  correct  positions,  for  both  the  zero- 
and  nonzero-offset  data.  Figure  4  shows  Gaussian  beam  migration  result  by  using 
the  same  zero-offset  data.  While  both  our  method  and  Gaussian  beam  migration  did 
well,  the  Gaussian  beam  migration  used  17  minutes  on  the  IBM  Risk  System/6000, 
and  our  method  used  7  minutes,  even  though  I  used  the  prestack  algorithm  for  the 
zero-offset  data. 

The  second  example  is  on  turning-wave  migration.  The  model  shown  in  Figure  5 
is  a  single  reflector,  with  a  segment  beyond  90  degrees.  The  velocity  function  is  the 
same  as  in  the  first  example,  shown  in  Figure  lb.  With  the  model  and  the  velocity 
function,  I  generated  a  data  set  with  an  offset  of  500  m,  shown  in  Figure  6a.  After 
migration,  the  reflector  is  imaged  to  its  correct  position  shown  in  Figure  6b,  even 
for  the  segment  beyound  90- degrees.  This  result  shows  that  my  method  can  handle 
lateral  variation  in  velocity  and  turning  waves. 

Smoothing  test 

The  velocity  model  shown  in  Figure  7  consists  of  six  constant-velocity  layers  and 
the  input  zero-offset  data  for  migration  shown  in  Figure  8  were  generated  with  this 
model  by  the  CSHOT  program (Docherty,  1987). 

Different  smoothed  velocity  functions  are  used  in  Gaussian-beam  migration  that 
is  implemented  by  Hale’s  program  (1992).  Figure  9  shows  the  data  migrated  with 
the  unsmoothed  velocity  in  Figure  7.  The  migration  result  is  poor:  low  resolution 
and  migration  artifacts  are  in  the  first  three  reflectors,  and  the  last  two  reflectors 
are  almost  invisible.  Figure  10  shows  the  data  migrated  with  a  smoothed  velocity 
(smoothing  parameters  ax  =  az  =  2.5).  While  it  is  much  better  than  Figure  9,  there 
are  still  some  artifacts.  Figure  11,  showing  the  data  migrated  with  a  more  smoothed 
velocity  (  smoothing  parameter  ax  =  az  =  5),  gives  an  improved  structural  image.  If 
too  large  smoothing  parameters  (ax  =  a2  =  20)  are  used,  the  migration  result,  shown 
in  Figure  12,  deteriorates.  One  can  see  that  the  corners  of  the  third  reflector  have 
become  too  sharp,  and  the  bottom  reflector  is  no  longer  fiat.  Figure  13  shows  the 
data  migrated  with  the  smoothed  velocity  by  the  windowed  averaging  method.  The 
window  size  is  selected  so  that  the  difference  between  the  smoothed  velocity  and  the 
original  one  is  about  the  same  as  that  in  Figure  11.  Both  the  damped  least-squares 
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and  the  windowed  averaging  method  improve  structural  images  through  smoothing 
the  velocity,  but  the  former  looks  a  little  better  than  the  latter.  There  is  incoherence 
on  the  vicinity  of  the  right  intersection  between  the  first  and  the  second  interface  in 
Figure  13  compared  to  Figure  11.  The  selected  parts  are  shown  in  Figure  14.  From 
this  result,  one  may  conclude  that  if  the  differences  between  the  smoothed  velocities 
and  the  original  one  are  the  same,  the  damped  least-squares  method  smooths  velocity 
a  little  more  effectively,  therefore,  gives  a  slightly  better  result. 

The  velocity  in  Figure  7  has  strong  lateral  variation.  For  this  unsmoothed  velocity, 
the  finite  difference  algorithm  for  the  eikonal  equation  6  will  break  down  because 
dr/dz  may  not  positive.  Therefore,  in  order  to  implement  the  Kirchhoff  migration 
of  this  paper,  the  smoothed  velocity  is  required.  The  least  smoothing  parameters  for 
this  requirement  are  ax  =  a*  =  14.  Since  the  velocity  is  slightly  oversmoothed,  the 
bottom  reflector  is  not  very  flat.  Furthermore,  like  other  Kirchhoff  migration  methods, 
the  discontinuities  in  the  bottom  event  of  Figure  8  make  strong  diffraction  smiles  in 
Figure  10  since  large  apertures  are  used.  In  contrast,  Gaussian  beam  migration  does 
not  suffer  this  problem. 

Figure  16  shows  the  differences  between  the  smoothed  velocities  and  the  original 
one,  which  is  calculated  by  using  formula  (20).  The  larger  the  smoothing  parameters, 
the  larger  the  difference.  This  difference  may  help  to  choose  suitable  smoothing 
parameters.  For  example,  the  difference  for  ax  =  ax  =  20  is  greater  than  0.1  from 
which  I  conclude  that  the  velocity  is  oversmoothed  and  then  reject  these  smoothing 
parameters. 


CONCLUSION 

The  Kirchhoff  integral  provides  a  powerful  tool  for  modeling  and  migration.  In  this 
paper,  a  finite-difference  algorithm  is  used  to  calculate  traveltimes  and  amplitudes. 
With  the  help  of  interpolation,  this  method  can  be  efficiently  implemented.  The  result 
on  the  synthetic  data  shows  that  this  method  can  handle  lateral  variation  in  velocity 
and  turning  waves.  One  limitation  is  that  the  Kirchhoff  method  cannot  efficiently  deal 
with  caustics  in  the  Green’s  functions  of  the  integration  operator  and  multiple-arrival 
times.  Another  is  that  this  method  needs  a  sufficiently  smooth  velocity  function. 
When  velocity  has  a  strong  lateral  variation,  the  eikonal  equation  cannot  be  solved 
by  this  method.  This  difficulty  may  be  solved  through  udng  smoothed  velocity. 
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APPENDIX  A:  GEOMETRIC-SPREADING  FACTOR  FOR  LINEAR 

VELOCITY 


Firstly,  I  assume  that  the  velocity  is 

v  (z)  =  r>0  +  oz. 

Following  from  formulas  in  Bleistein  (1986),  I  find 


vo 


x  —  x3  =  — : — r(cos  d  —  cos  </>), 
o  sin  d 


and 


<7  = 


From  Snell’s  law 


Vq  (cos  d  —  cos<Z>) 
a  sin2/? 

sin  jd  s\n4> 

vq  v(z )  ’ 


I  have 


d(f>  r(  z)  cos  d 
dd  co  cos  4> 

By  using  the  Snell’s  law  and  equation  (A-l),  I  obtain 

i>0  o  viz) 

X  —  Xg  =  —  cot  d - cot  <t>. 

a  a 


(A-l) 

(A-2) 


(A-3) 


(A-4) 
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Fig.  2.  (a)  Zero-offset  data  for  the  model  in  Figure  1.  (b)  The  Kirchhoff  migration 
of  the  data  in  (a),  with  velocity  model  in  Figure  lb. 
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Fig.  3.  (a)  Nonzero-offset  data  (offset=1000  m)  for  the  model  in  Figure  1.  (b)  The 
Kirchhoff  migration  of  the  data  in  (a),  with  velocity  model  in  Figure  lb. 
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FlG.  6.  (a)  Data  (offset=500  m)  for  the  model  in  Figure  5.  (b)  The  Kirchhoff 
migration  of  the  data  in  (a),  with  velocity  model  in  Figure  lb. 
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FlG.  8.  Zero-offset  synthetic  data  generated  with  the  velocity  model  in  Figure  7. 
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Fig.  10.  Gaussian  beam  migration  with  a  smoothed  velocity  model.  The  smoothing 
parameters  ax  =  az  =  2.5.  The  velocity  is  undersmoothed. 
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11.  Gaussian  beam  migration  with  a  smoothed  velocity  model.  The  smoothing 

parameters  ax  =  az  =  5. 


Fig.  12.  Gaussian  beam  migration  with  a  smoothed  velocity  model.  The  smoothing 
parameters  ax  =  az  =  20.  The  velocity  is  oversmoothed. 
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Fig.  13.  Gaussian  beam  migration  with  a  smoothed  velocity  model.  The  velocity  is 
smoothed  by  the  windowed  averaging  with  the  window  size  of  9  points. 
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Fig.  14.  Selection  form  Figure  11  and  Figure  13.  The  left  one  is  from  Figure  11  and 

the  right  one  is  from  Figure  13. 
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Fig.  15.  Kirchhoff  migration  with  a  smoothed  velocity  model.  The  smoothing 
parameters  ax  =  az  =  14.  The  velocity  is  slightly  oversmoothed. 
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FlG.  16.  The  relative  rms-differences  between  smoothed  velocities  and  the  true  ve¬ 
locity.  Four  solid  curves  denote  the  smoothed  velocities  by  the  damped  least-squares. 
The  dashed  curve  denotes  the  smoothed  velocity  by  the  windowed  averaging. 
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Taking  the  derivative  with  respect  to  3  and  using  equation  (A-3)  gives 

dx  __  vq  v(z)  d<t>  _  t’o  v2(z)  cos  3  vq(cos  3  —  cos  0) 

03  a  sin2  3  a  sin2  0  d3  a  sin2  3  av 0  sin2  4>  cos  0  a  sin2  3  cos  0 


Therefore, 


dx  _  a 
03  vq  cos  0  ’ 


(A-5) 


03  _  i’ocos0 
Ox  a 


From  equation  (4), 


03  03  Or  [dr]  1  03 

Oz  Ox  Ox  Oz  Ox tan  ^ 


The  above  formula  and  the  chain  rule  for  partial  differentiation  give 

03  03  Ox  03  Oz  03  r  cos {9  -  0) 

7TT  =  Tprrr  +  =  —  ( r  cos  6  +  r  sm  9  tan  0)  = - 

06  Ox  00  Oz  06  Ox  cos  0 

Finally,  substituting  equation  (A-6)  into  (A-7),  I  obtain 

i 


(A-6) 


(A-7) 


(A-8) 


When  the  velocity  is  a  linear  function  of  x  and  z,  there  is  a  rotational  transfor¬ 
mation  around  the  point  (xs,0)  such  that  the  velocity  is  a  linear  function  of  z  in  the 
new  coordinate  system.  Also  it  is  easy  to  show  that  equation  (A-8)  is  invariant  under 
the  rotational  transformation. 

APPENDIX  B:  SMOOTHING  IN  WAVENUMBER 

By  using  the  Euler’s  formula  from  the  calculus  of  variations,  we  change  (16)  into 
a  differential  equation  for  vs(x,  z): 


v,(x,  z)  -  v(x,  z)  -  a\  — ■  =  0. 


(B-l) 


Taking  the  Fourier  transforms  with  respect  to  x  and  2  gives  the  solution  in  the 
wavenumber  domain 

vs(kx,kz)  =  v(kx,kz)/{  1  +  ot2xk2). 

Similarly, 

vs(kz,kz)  =  v3{kx,kz)/(  1  +a2k2). 
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Vg  ( kx i  kz)  — 


v(kx,k2) 

(l+a2fc2)(l+a2fc2)- 


In  discretization,  the  second  derivative  can  be  approximated  by 

d2v3  ^  v3{x  4-  Ax,  z)  -  2i3,(x,  z)  4-  t\(x  -  Ax,  z) 
~d^  “  (Ax)2 

whose  Fourier  transform  with  respect  to  x  is 

4sin2(fcxAx/2)  _  s 
- - v^z^ 


(B-2) 


Using  this  formula  and  taking  the  Fourier  transforms  with  respect  to  x  and  z  in  (B-l) 
give  the  solution  in  the  wavenumber  domain 


Similarly, 


Vs{kx i  kz)  — 


vs(kz,kz)  = 


v{kx,  k2) 

1  4-  4(aI/Ax)2sini(fcrAx/2) 

_ va(kx,k2) _ 

1  4-  4(a*/Az)2  sin2(fcz Az/2)  ’ 


v  (h  k  )  = _ v[kz,kz) _  (B-31 

(1  4- 4(ar/Ax)2sin2(fcxAx/2))  (1  4- 4(ar/Az)2sin2(fczAz/2))’ 
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